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Abstract. In this paper, we will study the following geometric flow, obtained by Goldstein and 
Petrich while considering the evolution of a vortex patch in the plane under Euler's equations, 

Xt = -ksu ~ -k'^T, 
00 ^ 

C ^ ^ with s being the arc-length parameter and k the curvature. Perelman and Vega proved in 117] that 

C ^ this flow has a one-parameter family of regular solutions that develop a corner-shaped singularity at 

flnite time. We will give a method to reproduce numerically the evolution of those solutions, as well 
^ as the formation of the corner, showing several properties associated to them. 

(D 

Q 

in 

< 

Xt^-ksxi- -k^T, (1.1) 

S . ^ 

with s being the arc-length parameter and k the curvature, Tg = fcn. It was obtained 
I by Goldstein and Petrich in |llj : their motivation was the problem of the evolution 

^ of a vortex patch in the plane subject to Euler equations [TF . If the boundary is at 

CO least piecewise C^, the exact motion of the boundary of the patch satisfies 

S X.(a.,)^-^ri„( "°'"-^'°'"ll )x.(o-,,M»-. (L2, 

CN 

where a is the Lagrangian parameter, ujp is the vorticity and tq is an arbitrary pa- 
rameter whose choice does not affect the dynamics of the curve. 
L* Let us rewrite (1.2 1 using the arc-length parameter s; then 
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1. Introduction. In this paper, we will consider the following geometric flow of 
planar curves that can develop singularities at finite time 



^ Xt(s,t) = -// In ^^^^ X,(,s',i)ds', (1.3) 

with L(t) being the length at time t. We truncate this last integral by introducing a 
cutoff at s = ±A: 

X.(M)»-^ r-'l„( l|X(M)-X(„-,,)lh ^^(,,„,,, 

27r Js^A \ ro J 
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We expand X(s', t) and Xs(s', into powers of A = s — s', bearing in mind — T, 
-kT: 



Tg = fcn, ris 



1 



X(s') = X(s) + T(s)A + ^fc(s)n(s)A2 + • • • 
T(s') = T(s) + fc(s)n(s)A + ^[A;,(s)n(s) - A:2(s)T(s)] A^ + 



Introducing the Taylor expansions into ( 1 .4 1 , we can integrate term by term. If we 



represent the PDE for X as Xj = Un + WT, then, considering the leading terms of 
the expansions, we obtain the following approximations for U and W: 



U 
W 



'^kM (inA-i-ln(ro) 

^EA(ln A - 1 - In(ro)) + ^^k^s) Lk-\-\- In(ro) 

TT DTT \ 3 4 



We choose tq = Ae in order to have Wg 
approximation for X is 



kU [TD]. With that choice, the final 



,A3 



367r 



1 



(1.5) 



A Galilean transformation removes the term — ^T; then, the factor — is ab- 

27r ' ' 367r 



sorbed after a change of variable, getting (|l.ll 

Xt 



1 . 



A:,n- -fc^T. 
2 

Since we are considering planar curves, we can identify the plane where they live with 
C; denoting z = X and bearing in mind that Zss — ikzg, the last equation becomes 



1, 



(1.6) 



In this form, the local induction approximation preserves some of the basic conserved 
quantities of the exact vortex patch dynamics; for example, area, center of mass and 
angular momentum [5]. It also preserves Sq and, in particular, the length, which 
is not true in the vortex patch problem. Indeed, numerical calculations O |S] show 
that small bumps in the boundary of isolated vortex patches of constant curvature 
cause filamentation phenomena to occur, i.e., the ejection of thin filaments into the 
surrounding fluid. Nevertheless, for a small enough perturbation, the time at which 
filamentation appears can be made arbitrarily large and we can assume that the initial 
parametrization of the curve Sq,, and hence L, are time-independent. 

Flow (1.6 1 is time-reversible, because if z{s,t) is a solution, so is z{—s,—t). It 
is completely determined by its curvature, k{s,t), except for a rigid movement that 
changes with time and that can be fixed by the initial conditions. As shown by 
Goldstein and Petrich, k satisfies the modified Korteweg-De Vries (mKdV) equation 



-k^ks = 0. 
2 



(1.7) 



Numerical study of a flow that develops singularities at finite time 



3 



To relate the vortex patch evolution and the niKdV equation, Goldstein and Petrich 
followed previous ideas by Hasimoto [T^], who connected the nonlinear Schrodinger 
(NLS) equation with the motion of vortex filaments in , ideas which were extended 
by Lamb [H]. Later, Nakayama, Segur and Wadati 16J identified the connection 
between integrable evolution equations and the motion of curves in the plane and 
in M.^. More recently, Wexler and Dorsey [19] found that under a local induction 
approximation, the contour dynamics of the edge of a two-dimensional electron system 
can be described again by the mKdV equation. 

In [T7], Perelman and Vega proved the existence of a regular family of solutions 
for (|1.6| that develop corner-shaped singularities at finite time. Conversely, they also 



proved the existence of solutions of the mKdV equation (1.7 1 with initial conditions 
given by 



k{s, 0) = aS{s), 



a e 



(1.8) 



where d{s) is the Dirac delta function and \a\ is small enough. The corresponding 
initial condition for (1.6 1 is 



z{.s,0) 




(1.9) 



for some e+,9- e [0, 2tt), 6+ - 6 
1 looked i 

k{s,t) = 



2a. 



Perelman and Vega looked for self-similar solutions of (1.7) in the following form 

I ( s 



t > 0. 



(3i)i/3 \{3ty/3 
which leads to study the following ODE 

Uxx ^ xu + 2u^ = /-i, a: G M, /i € M, 
being /i an integration constant. We will only consider the case /i = 0, 

Uxx — xu + 2u^ = 0, a; G M, 
but the method developed in this paper can be easily implemented for fJ. =/= 0. 



Equivalently, the self-similar solutions for (1.6 1 are of the form 
z(s,i) =ti/^c^(^) , t>0, 

which leads to study 




s e 



(1.10) 



(1.11) 



(1.12) 



(1.13) 



(1.14) 



Bearing in mind all the previous arguments, Perelman and Vega proved the following 
theorems: 

Theorem 1.1. There is eo > such that if a? < eo, then there exist 9^ G [0,27r) 
and 10 an analytic solution of (1.141 such that if 



t > 0, 



(1.15) 



then 
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(i) z solves ( |1.6[ ) /or t > and 

z{s,t) - se'^*X[o,oo)(s) - se*" X(~oo,o](s) 

wii/i X being the characteristic Junction, 
(n) e+ -0- = 2a; 



< ct^l^; 



(Hi) The curvature k of z satisfies (1.10 I and (1.12). 

Theorem 1.2. There is eo > such that given any a, with < cq, there exists 



a bounded real analytic u{x,t) solution of (1.121 such that if 



k{s,t) 



2 / s 



t > 0, 



(1.16) 



then 



(i) k solves the mKdV equation (1.7); 

(•+00 

(a) I u{x)dx = a; 



(Hi) lim/c(-,t) = 2a5 in S' 



the space of tempered distributions. 
These theorems, which constitute the theoretical basis of this paper, guarantee the 



formation of a corner-shaped singularity at finite time for (1.6), provided that \a\ is 



small enough. Nonetheless, numerical simulations in subsection |3 .2| will give evidence 
that singularity formation happens also for any a G (— 7r/2, 7r/2), i.e., for parameter 
values outside the scope of Perelman and Vega's theory. 



The purpose of this work is to study the self- similar solutions of (1.7) from a 



numerical point of view, as well as the formation of their corresponding corner (1.9 1, 
going backwards in time from t — 1 until t — 0, because ( |1.7[) is also time reversible. 

Instead of developing a numerical method for (1.6) or (1.7), we rather consider 
the angle 



9{s, t) = e{-oo, t)+ k{s', t)ds', 



(1.17) 



and the PDE for the angle, obtained after integrating (1.7 1 once, 

1 



t{s,t) 



(1.18) 



Working with 9 has two main advantages: it allows to guarantee naturally \zs{s, t)\ 



for all t and we can preserve numerically the conserved quantity of (1.7) 
2 



1 p+co r+oo 

- u{s)ds= k{s,t)ds = 6+ - e-, 

^ J -oo J-oa 



by fixing 9{~oo, t) = 0^ and 0{+co, t) = 0+ , for all t. 

The structure of this paper is as follows: In section [2] we integrate (1.12), looking 
for admissible initial data {u{0) , Ux{0)) , such that the corresponding solutions satisfy 
u{x) —^ as X —^ oo; this must be carefully done, because those solutions are very 
unstable. 

The admissible pairs {u{0),Ux{0)) form a connected curve. Each point of this 
curve determines one solution for (|1.6p , ( 1 .7 1 and (1.18); hence, we have one-parameter 
families of regular solutions for (1.6), (1.7) and (1.18) that develop singularities at 
finite time. 
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In section |3] a spectral numerical method with integrating factor for (1.18 1 is 



developed, having truncated s e M to [saiSb], with <C —1, Sf, 3> 1. We impose 
the boundary conditions 9(sa,t) = 9^, 9(sb,t) = 9^, which is equivalent to fixing the 
tangent vectors of z at Sa and st- In subsection |3 . 1 1 we explain how to integrate 9{s, 1) 
in [sa, Sb], which involves an estimate of J^°° k{s, 1). Finally, numerical experiments 
are carried out in subsection 13.21 

In the exact problem, the energy k^{s,t)ds = cx) is a conserved quantity for 
all t; this infinite energy concentrates at s = as t — > 0, which causes z to develop 
the corner-shaped singularity. In our numerical experiments, the energy in [sci,Sb], 
'' k'^{s,t)ds, is finite, but, nevertheless, it keeps approximately constant and it also 
tends to concentrate at s = 0, as t — ^ 0. This fact shows that even after having 
truncated K to [sajS;,], the energy accumulation process continues to be stable. It is 
also remarkable the good accuracy with which we recover even 



for small t, hence approaching the Dirac delta function (1.8 1. The numerical results 
suggest that the accuracy of /c(0, t) could be improved arbitrarily by increasing the 
length of [sa, Sh], i.e., the energy of the system; it would be very interesting to prove 
analytically that we can recover the solution of the exact problem by making Sj, — s- — oo 
and Sb — > oo. 

In sectionjij we calculate the estimates of k{s, 1), as explained in subsection 



3.1 for a large set of admissible initial data (u(0), ^^^(O)) of (1.121, giving numerical 



evidence that 

/ u{x)dx ^ ("o^'o^j or, equivalently, / fc(s', i)fis' G (— tt, tt), Vt. 

In section [5] motivated by the original vortex patch problem, we append to the initial 
datum 9{s,l) a smooth function in such a way that the corresponding 2(s,l) is a 
closed regular curve without intersections. Numerical experiments with this new 
initial datum are carried out in subsection |5.1[ showing that the method developed 
in section [3] keeps z{s^t) closed for all t, preserving its inner area as well. It is also 
observed that closing z has no effect on the energy concentration process. 
Finally, conclusions are summarized in section |6] 

Besides the local induction approximation done by Goldstein and Petrich |llj . 
some other simplified models have been proposed to describe the vortex patch dynam- 
ics. In [6j, Constantin and Titi introduced a hierarchy of area-preserving nonlinear 
approximate equations, showing that the first of these equations, starting from arbi- 
trarily small neighborhoods of the circular vortex patch, blows up. Later on, Alinhac 
[T] considered a quadratic non area-preserving approximation for vortex patches with 
contours near the unity circle, obtaining an instability result at finite time. 

In the vortex patch problem, Chemin has proven in that, when considering 
a smooth initial contour, no finite-time singularities may happen (infinite length, 
corners or cusps, for instance), i.e., smooth contours stay smooth in time; later on, this 
result has been proven also by Bertozzi and Constantin . There is no contradiction 
between those results and our model, because Kenig, Ponce and Vega '13] have proven 
in M that finite energy solutions of the mKdV equation do not develop singularities; 
an equivalent of this has been proved by Bourgain 0] in the torus. In |17j . Perelman 
and Vega were instead considering infinite energy solutions of the mKdV equation, 
which caused the singularity to happen. 



Our experiments exhibit the energy concentration process in (1.6 1, both after 



having truncated z and after having closed it with a big loop. Since we consider finite 
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energy solutions, we cannot reproduce the corner-shaped singularity formation, but 
just approach it. It would be very interesting to simulate the evolution of this closed 
curve under the equations for the vortex patch. Although no singularities can be 
obtained, it could be still expected to see some energy concentration process taking 
place (obtaining "big curvatures" from "small" initial curvatures). 

Despite its good properties (conservation of area, center of mass, angular mo- 
mentum) and the fact that it does not create singularities, the main criticism of our 
model is that the arc-length parameter is preserved for a simple closed curve and, 
therefore, the total length of the curve remains constant. This does not happen in the 
vortex patch problem, where there are examples in which the length or curvature of 
the vortex patch boundary grow rapidly [11 16]. Therefore, the mKdV equation should 
be complemented by another equation for the evolution of s^. 



2. Integration of u{x). Let us come back to (1.12) 



= xu — 2u 



X e 



(2.1) 



In order to integrate this second-order ODE, we rewrite it as 



Vj. — XU — 2u^ 



(2.2) 



and use the fourth-order Runge-Kutta method, needing two initial data u(0) and 
UxiO) — ^{0). j{x) = u{x')dx' can be also integrate, without extra effort, by 
adding the equation 7x = u to the system (2.2 1. This idea will be used in subsection 



If u is a solution of (2.11, then, for a; ^ 0, u{x) has the same oscillatory behavior as 

0. Likewise, when x ^ oo, the solutions 



the solutions of the Airy equation u^ 



- xu ■ 



of (2.1) are characterized by being very sensitive with respect to small variations 
of the initial data, i. e., if u{x) oo, small changes in the initial data can make 
u{x) — *■ — oo and vice versa. Let us consider, for instance, u{0) — 0.024, as well as 
four possible choices for Ux{0) and let us integrate numerically ( |2.1[ ). 



u(x) for different initial data 




Figure 2.1. Dependance of u{x) on the initial conditions 



In figure 2.1 the solutions for Ux{0) = —0.018 and ^^,(0) = —0.017, plotted dis- 
continuously, tend respectively to — oo and oo, and the solutions for Ux{0) = —0.0175 
and Ux{0) = —0.0174, plotted continuously, tend respectively to — oo and oo, although 
we observe that for these last two values, the explosion happens a bit later. 

If we go farther with the process, between Ux{0) — —0.0175 and Ux{0) = —0.0174 
there exists a unique Ux{0) such that liiiix^coU^x) — 0. That value, obtained with 
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a bisection technique, is approximately Ux{0) — —0.0174881944. . . In figure 2.1 for 



scale reasons, we cannot distinguish clearly the oscillations in the real negative axis; 



in figure 2.2 having used the limit value for Mx(0), those oscillations are displayed. 



u(0) = 0.024, u (0) = -0.0174881944... 




Figure 2.2. u(x), with lim u{x) = 0. 



Because of the sensitivity with respect to initial data, we can calculate numerically 
only as many decimals of the correct Ux{0) as the machine precision allows us, which 
implies that we are only delaying the explosion time for u(x). We can say in an 
equivalent way that the solutions that make 



lim u{x) — 



(2.3) 



are highly unstable. Nonetheless, for our purposes, it is enough to consider u{x) = 

for a big enough x, because u{x) — > exponentially as a; — s- oo. Hence, since 

1 u{x)dx\ < cxD, due to the fact that the oscillations in x < cancel one another, 
condition ( |2.3[ ) is equivalent to | JJ*^ u{x)dx\ < oo. 




Figure 2.3. Pairs (u{0) , Ux {0)) , with lim u{x) = +oo (black) and lim u{x) = — oo (white) 



In figure |2T3| the points belonging to the region in black, when taken as initial data 
(it(0), Ux{0)) of (2.1 1, give us solutions for ( |2.1[ ) such that u{x) +oo, as a; ^ oo; for 
the points of the region in white, we have u(x) — > — oo instead. The boundary (figure 
2.4) between both regions are the admissible pairs (u(0),Ux{0)) whose corresponding 
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solutions satisfy (2.3 1. That boundary is a connected unbounded one-dimensional 
curve that divides the plane in two antisymmetrical halves. Each of the points of that 
curve determines one single self-similar k(s,t), one single self-similar 0{s,t) and one 
single self-similar z{s,t), which are solutions of (1.7l, (1.18) and (1.6l, respectively. 
Therefore, we have one-parameter families of solutions for (1.7l, (1.181 and (1.6 1 that 
develop a singularity at finite time. 




3. A numerical method for the evolution of (1.6 1. The k{s,t), 9{s,t) and 
z(s,t) we are considering are characterized by 

lim |z(s,i) - se'^^l = 0, lim |z(s, i) - se'*'" | 0. 

s — ^ — oo s — >+oo 

Taking 6{s,t) as the evolution variable, we can preserver naturally k{s,t)ds and 
ks(s, ^)l = 1- The PDE for 9 with boundary conditions is 

9t{s,t) = -9s,sis,t) - ^{esf{s,t), s e R, 

e{^oo,t)^e-, (3-1) 

(+C50,<) = 9+. 

Due to the difficulty of considering the whole M, we will restrict ourselves to s € [sq, Sb], 
with Sa <C —1 and Sh ^ 1 

9t{s, t) = t) - ]^{esf{s, t), S e [Sa, Sb], 

e{sa,t) = 9-, (3-2) 

e{sb,t) = e+. 

9{s,t) is not periodic in [sci,s;,], so we define 

~e{s,t) = e{s,t)~'^{e+ -e-)-e-, L = st-Sa, (3.3) 
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which is periodic and regular for alH > 0. (3.2 1 gets transformed into 



L 



S G [Sa,Sb], 



(3.4) 



We will calculate the numerical evolution of at + 1 equidistant points {sj} in 
[sa, Sb]- Denoting with some abuse of notation d{^,t) = {d{S^,t))^, we represent 



N/2-1 



(Sj - So) 



(3.5) 



where L ~ Sb — Sa, As — ^ and sj ~ Sa + J As. Therefore, (3.4 1 is transformed into 



1- \ 3 



L 



(e,i), (3.6) 



with ^ = — — 1. Since we are working with the frequency, the third derivative 
is transformed into a multiplier that can be absorbed by means of an integrating factor 



ei^,t) exp 



■ exp 



e+ - e- 



A 



The advantage of using an integrating factor is twofold: it allows to integrate exactly 
the linear part of (3.2 1, increasing the accuracy of the numerical results, and it relaxes 
considerably the time-step restrictions. 

We apply the fourth-order Runge-Kutta in time with integrating factor as de- 
scribed in [TH]. Denoting 



e°(s) = 0(s,i"), 
e"(s)«r(s)^0~(s,i"), 

we have 



Mo 

Mo 
e(^)(0 



e"(0 ~r(C) = ^(e,i"), r = t° + nAt, 



exp 



= exp 




A 



(0 



(6, 



e(-^)(0 



L 



(6, 



"(^)e(-)(0 



e+ - e- 



(0, 



10 



Francisco do la Hoz 




The symbols A and V denote respectively the direct and inverse fast Fourier transforms 
(FFT) 0. 

Finally, we force in every time step 



@n+i ^ (e«+i)v ^ sR((e"+i)V) 
e"+i(so) - e"+HsAr) = 0. 



(3.7) 



Rounding 0"+"'^ to zero at the boundary points Sq = sq and si, — sn avoids the 
accumulation of little errors of order 0(10^^). 



by 



3.1. Computation of 9{s, 1). Bearing in mind (1.10 1, at t — I, 9{s, 1) is given 



(»,1)- 9(0,1)= / <■(»',!)*' = , 



where 



7(2^) — / u{x')dx'; 



(3.8) 



(3.9) 



hence, = u and (2.2) can be generalized to 




(3.10) 

with initial data m(0), Ux{0) = v{0) and 7(0) = 0. (m(0),u(0)) must be an admissible 
initial pair for (2.2 1, i.e., such that the corresponding u{x) satisfies lima;^oo m(x) = 0. 

From now on, wc wi ll cho ose m(0) = 0.72, v{0) = 1.1601860809647328. In figure 
3?T| we have integrated ( |3.10| in a; e [-80,20], with |Aj:| = lO"'^. For x > 9.204, we 
have taken k{x) = 0, because \u{x)\ < 10^^^. We obtain 9{s,l) from (3.8l; hence, 
s — 3^^'^x and s £ [sa,Sb] — [—115.38,28.84]. Since we have one degree of freedom, 
we fix 9+ = 9{sb,l) = 0. 

As we observe in the lower part of figure 3.1 0(s, 1) tends to 0~ extremely 
slowly as s — > —00. Nevertheless, we can determine lims^_oo ^(s, 1) with high ac- 
curacy, as the mean value of the first maximum and the first minimum of 9(s,l), 
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for s > Sa- In our example, the first maximum takes place at s = —115.31, be- 
ing 6'(— 115.31, 1) = —2.8403; the first minimum takes place at s = —114.80, being 
6'(-114.79,l) = -3.0191. Thus, 

/ + 00 
k{s', l)ds' = 9+ -0- ^ 2.9294, 
-oo 

quantity plotted in the lower part of figure |3.1| with a thinner stroke. In section |4j 
we will approximate the value of J^°^ k{s' , l)ds' as a function of the admissible pairs 

Ko),u'(o)). 



e(s, 1) 




-120 -100 -80 -60 -40 -20 

s 

zoom of 9(s, 1) 




-2.9 - 



-116 -114 -112 -110 -108 -106 -104 -102 -100 



Figure 3.1. 9(s,l). In the zoomed image, we have plotted with a thinner stroke the mean 
between the first maximum and the first minimum of6{s, 1), for s > Sa, which gives a good estimate 

ofe-. 



The 9{s,l) obtained after integrating (3.101 is not directly suitable as an initial 
datum of (3.2), because 6{sa, 1) 0^ and (eq. 3.3l 

0{s,l)^0{s,l)-'-^{0+ -0-)-e-, i = ,5fa-S„ (3.11) 

is not periodic in [sa, Sb\- Hence, we have to modify 0(s, 1) as follows: 

• We chose A^i G N and divide [sa, Sb] into A^i + 1 equidistant points {sj}, in 
such a way that G {sj}- 

Sj^Sa+jAs, As=*^--^, j = 0,...,iVi. 

iVi 

We have taken Ni = 2000; hence As « 0.07211 and Sieoi = 0. 

• We find the Sj with the lowest index, such that Sj > Smin, where Smin is the 
first minimum point of 0{s,l) in (sa, oo), and Sj satisfies that 0{sj,l) > 0~ 
and k{sj, 1) = 0s{sj, 1) > 0. We name it Sjoint- Here, Sjoint — si3 ~ —114.51, 
^(sjomt, 1) « -2.9150 and fc(s,o,„t, 1) « 0.5479. 

• We redefine 0{s,l), appending at s = Sjdnt £^ smooth function that has a 
first order contact with 6{s, 1) and tends exponentially to 0^ for s < Sa- The 
final expression for 0(s, 1) is 
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joint 7 



1) - 0-] exp 



joint 1 ^joint^ 



-5 joint J 



S ^ Syo, 



joint J 



where, with some notational abuse, 0(s,l) stands for the original and the 
corrected functions. 

We evaluate the new 9{s, 1) at the former iVi + 1 points {sj}, as well as in 
some new equidistant points {sj} outside [sa,S6]- It is important that the 
final number of {sj} is + 1 = 2" + 1, for some n e N, in order to apply 
FFT efficiently. In the new set {sj}^, we define — sq, — 



The new function 9{s,l) satisfies 9{sa) — , 9{sij) — 6^ and makes (3.3 1 periodic. 
We improve its regularity by applying a smooth spectral filter to (3.3 1 



0(e,l) = e(e,l)exp 



-10 



2.5ICI 

N 



25 



(3.12) 



where e{^, 1) = I))' 



Corrected e(s, 1 ) 



——mmmmmmmM 

-200 -150 -100 -50 50 100 

s 

Zoom of 9(s, 1 ) and corrected 0(s, 1 ) 




Figure 3.2. In the upper part, the final corrected 9{s, 1) is plotted. In the lower part, both the 
original and the corrected 8{s, 1) are plotted and zoomed near the joint (marked with a circle); the 
corrected 6{s, 1) appears with a thicker stroke. 



In figure [3T2| we have taken N ~ 4096. In the upper part, we have plotted the final 
corrected 9{s, 1), with s G [—226.14,69.23]. We observe two long constant segments 
at the extremes of 9{s, 1), because of the rather large choice of N; this is convenient 
to avoid periodicity phenomena, since the exact 9{s, 1) is not periodic. In the lower 
zoomed part, we have plotted both the original and the corrected 9{s, 1), highlighting 
the point s = Sjoint with a circle. 



In its new form, it is immediate to obtain 9s{s,t) from 
derivation. Using again representation (3.3), we get 



9{.3,t) 



S - Sg 

L 



-9-)-. 



L 



{s,t), through spectral 



L^Sb-Sa- (3.13) 
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Although the process to obtain 6{s, 1) may look rather artificial, if we calculate 9s{s, 1) 
through (3.131 and compare it with the k{s, 1) obtained from integration of (3.101 and 



(1.10 1, we get the following error table: 



max |6's(s,l) - fc(s, 1)1 < 10" 

se[-113.3,oo) 

max |6ls(s,l) - fc(s, 1)1 < 10" 

sG[-110.6,oo) 

max \eJs,l)-k(s,l)\<W 

se[-108.6,oo) 



Thus, fc(s, 1) is recovered from the corrected 9{s, 1) with high accuracy, except near 
the joint, Sjomt — —114.51. This is graphically illustrated in figure 3.3 observe that 
the support of Osis, 1) is now finite. 



G (s, 1) 




-200 -150 -100 -50 

s 

Zoom of 8 (s, 1 ) and k(s, 1 ) 




Figure 3.3. In the upper part, the derivative of the final corrected 9{s, 1) is plotted; observe 
that its support is finite. In the lower part, the derivative of the corrected 9{s, 1), plotted with a 
thicker stroke, is compared with the original k{s, 1). 

z(s, t) can also be immediately recovered from 6(s, t), except for a rigid movement, 
because Zs{s,t) — exp{i9{s,t)). To fully determine z{s,t), it is not complicate to see 
that 



z{0, t) = -2{-itf'^ [m'(0) + u2(o)] z^(0, t), 
where Zg{Q,t) is a time-independent constant with unit modulus. 



(3.14) 



3.2. Numerical experiments. The aim of the method we have developed is to 
try to describe numerically the formation of the singularity in the self-similar solutions 



of equation (1.6 1 satisfying (1.10 1. This corner-shaped singularity happens at finite. 



time; indeed, at < = 0, we have 

z{s,Q) 



'zQ + se^^*, s>0, 
Zo-Hse*"", s < 0. 



Equivalently, when t —t Q, the curvature tends to a Dirac delta function 



k{s, 0) — aS{s), a 



+ 00 



u{x)dx. 



(3.15) 



(3.16) 
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This happens because these self-similar solutions have infinite energy 

f + C30 



k (s, t)ds = oo, 



and it tends to concentrate at s = as i 0. 

In our numerical experiments, we are not considering the whole M, but s £ [sa, Sf,]. 
At s = Sa and s = Sb, we fix the tangent vector of z, i.e., 9{sa,t) = 9~ , 6{si,,t) = 0~^, 
for all t. As we observe in the upper part of figure |3.3[ the support of the initial 
k(s,l) of our numerical experiments is finite; hence, our numerical solutions have 
finite energy at i = 1: 

+ 00 /■Si 

k'^{s,l)ds^ / fc^(s, l)(is < oo. 



We have executed the method with the initial 9{s, 1) plotted in figure 3.2 



m(0) = 0.72, u'(0) = 1.1601860809647328, 

u{x) integrated when x € [—80,9.203], 

s e [-226.14, 69.23], N = 4096, As = 0.07211. 



(3.17) 



In order to measure the quality of our results, we analyze two quantities: The evolution 
of the energy at s € [sa, si,] 



k^{s,t)ds, ie[0,l], 



and the curvature at the origin s = 0, k{0,t). From (1.10 1 



(3.18) 



(3.19) 



(3i)i/3' 

When considering the non-truncated problem with s G M, since all the infinite energy 



tends to concentrate at s = 0, the amount of energy (3.18 1 in s G [sa, Sb] grows up as we 



Jk^(s, t)dt, At = -10" 




0.2 



Jk^(s, t)dt, At = -5- 10 



0.2 



0.4 0.6 



Figure 3.4. Numerical evolution of the energy in s £ [sa, with initial data from ( |3.17| . Un- 
like in the non-truncated case with s £ M, the energy in the interval s S [sajSj] keeps approximately 
constant . 

approach t = 0, tending to infinity. On the contrary, in our numerical experiments. 
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we have observed that the energy in the interval s G [sa,Sh] keeps approximately 
constant after having fixed 6{sa,t) — 6^ and 6'(sb,t) = 9^; hence, we are preventing 
the energy outside the interval s € [sa, Sb] from entering. 

In figure 3.4 we have plotted the energy (3.18) in s G [sa, Sb] for At = —10^"* and 
At — —5 • lO"'^; remark that the energy conservation improves by diminishing \ At\. 

On the other hand, the accuracy of fc(0,t) is rather poor for small t (see figure 
3.5). Numerical experiments show that the way to recover the curvature k(0,t) with 





Figure 3.5. Numerical curvature at s = 0, with initial data from ( |3.17| and At = —5 ■ 10~^; 
the right-hand side figure is a magnification of the left-hand figure; the exact value ( |3.19| is plotted 
with a dotted stroke. 



bigger accuracy near t = is not to diminish \At\ or As, but rather to introduce more 
energy to the system, i.e., to lengthen the support of the initial datum fc(s, 1). To 
illustrate it, we have executed the method also with the following initial data: 



u(0) = 0.72, u'(0) = 1.1601860809647328, 

u{x) integrated when x in [—400,9.203], 

s G [-872.27, 309.22], N 16384, As = 0.07211 



(3.20) 



and 



u(0) = 0.72, u'(0) = 1.1601860809647328, 

u{x) integrated when x in [—800,9.203], 

s G [-1744.62, 618.36], N = 32768, As = 0.07211. 



(3.21) 



Thus, we have made the initial support of k{s, 1) approximately five and ten times as 
big respectively, but with the same As as in the previous experiment. 

As we see in figure [3^ we have been able to recover k{0,t) until much smaller t. 
Remark that the results are better in the case where the support of fc(s, 1) is bigger 
(lower right-hand side box). On the other hand, the observations about the energy 
are valid again (see figure 3.7). 

In figure 3.8 we show the evolution of fc(s, t) in function of space and time, for the 
simulation with initial data ( 3.20| ). Clearly, the support of fc(s, t) tends to concentrate 
into s = with linear velocity, as t — s- 0. 

Once we have obtained 9{s,t), it is immediate to recover z{s,t), from 9{s,t), 
making Zs{s,t) = exp{i9{s,t)) , except for a rigid movement determined by ( 3.14[ ). 
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Figure 3.6. Curvature in s = 0. The upper part graphics correspond to the experiment with data 
\3.20\ ; the lower part graphics correspond to the experiment with data { 3.21\. T he right-hand side graphics 
show t £ [0, 0.01], i.e., time values close to the singularity; the exact value l |3.19| is plotted with dotted line. 



! k^(s, t)dt, A t = -10 



J k {s, t)dt, A t - -5-10 ■ 






Figure 3.7. Numerical evolution of the energy. The upper part graphics correspond to the 
experiment with data l|3.20| . the lower part graphics correspond to the experiment with data ||3.21[l. 





Figure 3.8. Numerical evolution of k{s,t). The support of k{s,t) tends to concentrate into 
s = with lineal velocity, as t — > 0. 



Nevertheless, from a numerical point of view, it is better for small times to fix z(sb, t) = 
z(sh, 0). In figure 3.9 we have superimposed the graphs of z(s, t) in a neighborhood of 
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s = 0, with initial data (3.201, at times t=l,t = 0.5, t = 0.1, t = 0.01 and t = 0.001, 



having fixed Zs(0,i) — 1. We can clearly appreciate the self-similar character of the 
solutions. 



z(s, t), with t IE {1 , 0.5, 0.1 , 0.01, 0.001} 




Figure 3.9. Numerical evolution of z{s,t) for different times, with initial data |3.20| l. The 
self-similarity is patent. The thick dark curve corresponds to t = 0.001. 



In figure |3.10| we have plotted with a thicker stroke the admissible pairs for 
which z{s^t) has no self-intersections, highlighting with a circle the pair that we have 
used in our experiments. 



Figure |3.9| and figure 3.10 explain why we have chosen 
(m(0),u'(0)) = (0.72,1.1601860809647328). Indeed, we have taken an admissible pair 
for which z{s,t) has no self-intersections, but such that it is near the situation when 
self-intersections happen; because of that, J k{s,t)ds > 2.92, i.e., a value near tt. 
Theorems 1 1 . 1 1 and 1 1 . 2| only guarantee the formation of the singular itv fo r small values 

give evidence 



of k{s' ,t)ds' — 6^ — 6 , but the previous results and figure 
that they also true for any value of 6^ — 6^ G (— tt tt — £),0<e^l (cf. section 
[2]). In any case, all admissible data on the dark, 's' shaped curve in figure 3.10 give 
self-similar corner solutions. 



3.10 




-2.5 -2 -1.5 



Figure 3.10. The thick dark curve shows the admissible pairs (m(0), Ma; (0)), such that z{s,t) 
has no self-intersections. The circle indicates the pair that we have considered. 



Before finishing this section, let us underline four main conclusions: 
• The finite energy in [sajSb] keeps approximately constant after fixing the 
tangent vector of z{s,t) at s = Sq and s = Sf,. 
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The support of k(s,t), and hence, the energy, tend to concentrate into s — 
with lineal velocity, as i 0. 

k{0, t) is well recovered, even for small t. To recover fc(0, t) for smaller times, 
it is necessary to consider a bigger support of fc(s, 1). 

Our numerical results generalize theorems |1.1| and |1.1| because they give 
evidence that the formation of the corner singularity happens for any value 
of 61+ -61- e (-7r + £:,7r-£), <£< 1. 



Hence, we are approximating a Dirac delta function ( 1.8 1, but we will never been able 



to create a singularity, because the energy is finite. In the non-truncated problem, 
with s e M, the energy is infinite for all t and it tends to concentrate at s = as 
t 0, which causes the singularity (1.8 1. 



Because of the parallelisms between the truncated and the non-truncated problem, 
we can say that our numerical solutions reproduce the non-truncated problem from a 
qualitative point of view, suggesting that we could approach the Dirac delta function 
as much as we want, by increasing the initial support of fc(s,l). It would be very 
interesting to prove analytically that we can recover the exact problem by making 

Sa "> — CXD, Sb — > +00. 



4. j k{s' ,t)dt, in terms of u{0) and Ux{0)- As seen in section 3.1 obtaining 
9{s, 1) for our experiments involves an estimate of k{s' , l)ds'. This integrate is 
hence determined by one admissible pair of initial data {u{0),Ux{0)) of (eq. 2.11, i.e., 
one point of the curve in figure |2.4[ 

We have integrated (3.101 in x e [—80,20], \Ax\ — 10^^, for a large set of admis- 



sible pair s (w(0), Ua;(0)) taken as initial conditions. Following exactly the same steps 
as in 3.1 we have estimated k{s' , l)ds' . 

In figure 4.1 each admissible pair (u(0),Ux(0)) appears together with its corre- 
sponding estimate of k{s' , l)ds'; the curve we have plotted in this way is anti- 



4 ^ 

3 

2^ 




u(0) 



Figure 4.1. Integral of the curvature. The left-hand side shows the admissible pairs 
{u{0),Ux{0)), together with their corresponding integral fl^^k{s',l)ds'. The right-hand side is 
a side view of the left-hand side. 

symmetric with respect to the origin. 

For the set of admissible data considered, our numerical estimates satisfy 

/ + 00 
fc(s', l,u{0),u'{0))ds' < vr + 5.6 • 10"^ (4.1) 
-OO 



Numerical study of a flow that develops singularities at finite time 



19 



Hence, there is strong numerical evidence to accept that 



+00 



u{x)dx e 



TT TT 

2' 2 



4- 00 



fc(s',t)ds' e (-7r,7r), \/t. 



(4.2) 



For values of the integral in [—3 — e,3 + e], < £ ^ 1, it is evident from the right- 
hand side part, which represents the first and third components of the curve on the 
left-hand side, that there is only one corresponding initial condition for (2.1). The 



values of the integral seem to converge exponentially to ±7r. We are prone to think 
that there is one single initial condition for each value in (— tt, tt), although if we have 
exponential convergence, it will be much more difficult to give numerical evidence for 
the values nearest to ±7r. 

5. Using a simple closed z(s, 1) as initial data. Coming back to our second 



experiment, with initial data (3.201, and motivated by the vortex patch problem, it 



is interesting to see what happens if we add a big loop to z(s, 1), so that it becomes 
a simple closed curve without intersections. For that purpose, we have appended 
a smooth function to 0(s, 1) at s = Sa, such that the new value for 6^ — 9^ is 
0^ — 9^ = 2ti. This function is a rescaling of 



Jo 



se [0,1], 



(5.1) 



with 




-f3 



(x — a)(l — x) 



[0,a)U{l}, 
s e (a,l). 



(5.2) 



being both ip and ^' are regular in [0, 1]. 

In what follows, we have rotated z in such a way that Zs{sa) — 1. If we plot the z 




Figure 5.1. 2(5,1) corresponding to the initial data | |3.20[ l before being closed 



corresponding to initial data (3.201, we observe in figure 5.1 that the lower asymptotic 
line of z is much longer than the upper one. That is why we have defined '(/'(s) being 
?/;(s) = in a certain [0,a], in order to prolong the upper asymptotic line of z before 
closing z by means of a loop; the choice of a determines the additional length that 
we add to the upper asymptotic line, whereas the choice of (3 determines the length 
of the loop. We have to adjust both parameters wisely, so that z becomes smoothly 
closed. 



20 



Francisco do la Hoz 



Remembering that 0{s, 1) was defined in such a way that 0(sb) — 0, the new 
expression for 0{s,l) looks like this: 



is,l)={27r-9{s,,l) f^^^\ (5.3) 



(s, 1), S e [Sq,S;,] 

Vl/(1) ^3(5,-5,)^ 

i.e., we have made four times as big the length of the new [s^, Sb] and our new si, is 
Sb = Sb + 3(s6 — Sa); the new amount of nodes sj is AN + 1 = 65537, so that As 
remains unchanged. Notice that, with some notational abuse, 9{s, 1) refers to both 
the old and the prolonged functions and Sb to the old and the new right extreme. The 
new 9 satisfies 9{sb, 1) — 9{sa, 1) = 27r. This is a necessary but not sufficient condition 
in order that z may be closed. To determine a and (3, we proceed as follows: 

If a is fixed, by a bisection technique, we determine a f3 such that ^(^(sb, 1)) — 0. 
Then, for some values of a, 3i{z{sb), 1) will be positive, while for some other values, 
it will be negative. For instance, for two choices of a, a = 0.1 and a = 0.2, we have 
applied the bisection technique to find the corresponding P in [0.5, 20]. At the end of 
the process, z{sb, 1) is not exactly real, but its imaginary part is negligible 



a = 0.1 



a = 0.2 



P = 5.492976875712412 
^2(56,1) « 209 + 5i • 10-12 

= 2.824543846424351 
z{sb,l) ~ -146-M0-12. 



Now, we have to apply the bisection technique to a. For example, for a — 0.15, we 
have 13 = 3.970063345682464 and z(s(,,l) « 32 - 2i • lO'^^ That means that the 



Prolonged e(s, 1 ) 



1000 2000 3000 4000 

s 

Original and ciosed z(s, 1) 




Figure 5.2. New 9(s, 1) and z{s, 1) 
correct a G [0.15,0.2]. After several iterations, the double-bisection algorithm yields 



(a,/3) = (0.158902181218767,3.741222383167766) 
z(sfc,l) « {2 + 3i) ■ 10-12. 



(5.4) 



Since s G [—872.27,3853.69] and z{s,l) is parameterized respect to the arc-length 
parameter, the length of the curve is 4725.96. In figure 5.2 we can see the prolonged 
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version of 6{s,l), as well as z(s, 1), with and without the loop. We have closed z{s,l) 
almost in a perfect way, because z{sh, 1) — z{sa, 1) = 0(10"^^). 



When z is closed, Zs(0, t) is obviously no longer constant and (3.141 does not hold 



any more; in fact, z{s, t) will rotate, so we need to give the evolution of a point and its 
angle for all t, which is quite straightforward. We just have to integrate the following 
two ode's in t for some sq: 



zt{so,t) = -Zsss{so, t) + -Zs{so, i)zss(so, t). 

The right-hand side of the first equation is known, so we can update 9{so, Then, 
bearing in mind that Zs(so,i") = exp(i0(so, the right-hand side of the second 
equation is also known and we update z(so,t"). 

In any case, since this paper aims at illustrating the energy concentration process, 
we will not rotate z, in order to compare the situations with and without the loop. 

log,Jz(s. , t) - z(s . t)|, Ate{-10-*, -510-=} 




Figure 5.3. Error when closing z{s,t). The error is smaller for At = —5 • 10 ^. 

5.1. Numerical experiments. We have executed the method with the new 
9(s,l). Remark that the evolution variable is 9, not z, so the method does not 
guarantee a priori that z{sa, t) = z{sb, t), for all t. Nevertheless, the curve keeps closed 



closed z(s, t) 



zoom of closed z{s, t) 





Figure 5.4. z{s,t) and zoom of z(s,t), for t £ {1,0.5,0.1,0.01,0.001} . The thick dark curve at 
the right-hand side corresponds to t = 0.001. 



with great accuracy, as we can see in figure 



5.3 



where we have plotted the decimal 
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logarithm of \z{sb,t) — z{sa,t)\, for At small enough. Indeed, for At = —5 • 10^^, 
k('Saj^) ~ z(sii,t)\ < 10^^ until very small times and, since the length of the curve is 
approximately 4726, that means a relative error smaller than 10~^. 

On the other hand, the loop does not have any visible influence on the formation 



of the corner-shaped singularity for z(s,t). In figure 5.4 we have plotted z(s,t), with 
Zs(0,t) = 1, for t = 1, 0.5, 0.1, 0.01 and 0.001. In the left-hand side part, all the 
plots seem to be a single one; nevertheless, after a potent zoom, the formation of the 
singularity at s = is clearly exhibited. In fact, the right-hand side of the figure is 
not ocularly distinguishable from figure [3^9] 

From the left-hand side of figure |5.4[ it is obvious that the enclosed area remains 
almost constant, for all t. We can check this by means of the well known formula: 



Area : 



{xdy - ydy) 



[x{s)y'{s) - x{s)y{s)]ds, 



(5.5) 



which is a direct application of Green's theorem. Figure [575| shows that the enclosed 
area is preserved with great accuracy. 



Area enclosed by z 




Figure 5.5. Preservation of the area enclosed by z 



6. Conclusions. In this paper, we have consider a geometric planar fiow (eq. |l.l[ ) 
found by Goldstein and Petrich |llj . while considering the evolution of a vortex patch. 
Perelman and Vega proved in |17j that it has a one-parameter family of solutions that 
develop corner-shaped singularities at finite time. We have studied those solutions 
from a numerical point of view, trying to reproduce the singularity formation. 

The flow can be characterized by the equation of its curvature, fc, which is the 
modified Korteweg-de Vries (1.7 1 equation. Nevertheless, we have found the angle 



= J kds' to be the most adequate evolution variable, because it allows to preserve 
naturally the arc-length parameter and the angle formed by the two asymptotic lines 
of z. Due to the difliculty of considering the whole M, we have taken s e [sa, Sb]- The 
main result is that, even after fixing 6(sa,t) = 9~ , 6{sb,t) = 0~^, i.e., after fixing the 
tangent vectors of z at Sa and Sb, we are still able to approximate the formation of 
the corner-shaped singularity. Indeed, the energy J^'' k'^{s',t)ds' < oo keeps approxi- 
mately constant and it tends to concentrate at s = 0; moreover, we recover with great 
accuracy fc(0,t), even for small t, tending /c to a Dirac delta function. The numerical 
results suggest that the accuracy of fc(0, t) could be improved arbitrarily by increasing 
the length of [sa,st,]; it would be very interesting to prove analytically that we can 
recover the solution of the exact problem by making — oo and s;, — > oo. 

The process of obtaining the initial 6 involves an estimate of the conserved quan- 



tity Jj;^ k{s',t)ds'. We have given numerical evidence that | k{s' ,t)ds'\ < tt 



+00 
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Motivated by the vortex patch problem, we have also considered a regular simple 
closed curve as initial datum, by appending a regular function to 0(s, 1) in such a 
way that z{s,l) is closed by a big loop. The method preserves the area enclosed 
by the curve with great accuracy and, more interestingly, the loop appears to have 
no influence on the energy concentration process; nevertheless, no singularities can 
be expected because the energy of our curves is finite [13l H] . This is in agreement 
with the vortex patch theory, where it is known that no singularities may arise from 
smooth initial contours [SJIS]- In [IZ], Perelman and Vega considered instead infinite 
energy solutions of the mKdV equation, which causes the singularity to happen. 

One of the shortcomings of our model is that the arc-length parameter is preserved 
for a simple closed curve and, hence, the total length of the curve remains constant. 
In the vortex patch problem, on the contrary, there are examples in which the length 
or curvature of the vortex patch boundary grow rapidly [Ull]- Therefore, the mKdV 
equation should be complemented by another equation for the evolution of Sq. 

7. Acknowledgements. The author would want to express his gratitude to 
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